%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Solve for Optimal Contract %%%
%Applies in general case with maturity delta>0
%Baseline is obtained via delta=0
%We solve the model given q and will insert the optimal q determined in
%OptContract.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Solution Procedure as in OptContract.m

%Solve model without screening moral hazard to obtain solution at the
%boundary V^B
SB;

%Boundary conditions for F' at V^B

if (aB<0.000001)

        dFStart=(FB-phi*(gamma-r))/phi;
        

end
    
right=1;
    
if (kappa*q<VB)
        
        right=0;
        
end

%%%%%%%%%% Solve HJB %%%%%%%%%%%%%%%
sol = ode15s(@ffunODE,[VB max(kappa*q,VB+0.001)+0.001*(kappa*q==VB)+0.001],[FB dFStart+tolerance-2*tolerance*right],options, r, gamma ,Lambda, phi,q,aBar,delta);


% Caclulate a, W, c, and F'' for all points in the state space
for i=1:length(sol.x)
    
    V=sol.x(i);
    F=sol.y(1,i);
    dF=sol.y(2,i);
    a=(F-dF*(sol.x(i)+phi)-phi*(gamma-r))/phi;
    a=min(F/phi,a);
    
         if (a<0)
            
            a=0;
          
            
        end
        
        if (a>aBar)
            
            a=aBar;
           
            
        end
        
        lambda=Lambda-a-q;
        W=phi*a;
        dV=(gamma+delta+lambda)*V-W;
        dF2=-(gamma+delta-r)*dF/dV;
        dW=-dF2*(V+phi)*dV;
        c=(gamma+lambda)*W+0.5*phi*a^2-dW;


        if(a>=F/phi-0.001)

            dW=dF;
            c=1;
        
        end
        
        vA(i)=a;
        vW(i)=W;
        vC(i)=c;
        vV(i)=V;
        vdV(i)=dV;
        vdF2(i)=dF2;
        
        
        
end

%Calculate time-0 values
[m,ind]=min(abs(sol.x-kappa*q));
a0=vA(ind);
aBOut=vA(1);
c0=vC(ind);
F0=sol.y(1,ind);
dF0=sol.y(2,ind);


%Calculate value of beta at V^B
betaBoundary=(betaB+0.00001)+(vW(1)-WB)/FB;

%Calculate Price of loan L(V), the retained share beta(V), and the expected
%time to default tau
solQuant = ode15s(@ffunQuant,[VB max(VB,kappa*q)+0.001],[LB betaBoundary tauB],options,mu, r, gamma ,Lambda, phi,q,sol.x,vA, vC,delta,VB,face,tauB,vA*phi,sol.y(1,:));

%%%%%%% Create equally-spaced grid of state space 
%%%%%%%%% Evaluate Quantities on this grid
quantGrid=linspace(VB, max(VB,kappa*q)+0.0001,10000);
vQuant=deval(solQuant,quantGrid);
quantGrid=linspace(VB, min(max(VB,kappa*q)+0.0001, max(sol.x)),10000);
solEval=deval(sol,quantGrid);


[m,ind]=min(abs(solQuant.x-kappa*q));
D0=solQuant.y(1,ind);
beta0=solQuant.y(2,ind);
tau0=solQuant.y(3,ind);

vD=solQuant.y(1,:)';
vbeta=solQuant.y(2,:)';
vbeta(1)=betaB;
vtau=solQuant.y(3,:)';

vbetaFine=vQuant(2,:);
vbetaFine(1)=betaB;